Sick-Sicker case-study  Supplementary Material to: A need for change! A coding framework for improving transparency in decision modeling

The Sick-Sicker model

In this case-study, we perform a cost-effectiveness analysis (CEA) using a previously published 4-state model called the Sick-Sicker model (E. A. Enns et al. 2015). In the Sick-Sicker model, a hypothetical disease affects individuals with an average age of 25 years and results in increased mortality, increased treatment costs and reduced quality of life (QoL). We simulate this hypothetical cohort of 25-year-old individuals over a lifetime (i.e., reaching an age of 100 years old) using 75 annual cycles, represented with n.t. The cohort starts in the “Healthy” health state (denoted “H”). Healthy individuals are at risk of developing the illness, at which point they would transition to the first stage of the disease (the “Sick” health state, denoted “S1”). Sick individuals are at risk of further progressing to a more severe stage (the “Sicker” health state, denoted “S2”), which is a constant probability in this case-study. There is a chance that individuals in the Sick state eventually recover and return back to the Healthy state. However, once an individual reaches the Sicker state, they cannot recover; that is, the probability of transitioning to the Sick or Healthy states from the Sicker state is zero. Individuals in the Healthy state face background mortality that is age-specific (i.e., time-dependent). Sick and Sicker individuals face an increased mortality expressed as a hazard rate ratio (HR) of 3 and 10, respectively, on the background mortality rate. Sick and Sicker individuals also experience increased health care costs and reduced QoL compared to healthy individuals. Once simulated individuals die, they transition to the “Dead” state (denoted “D”), where they remain. Figure shows the state-transition diagram of the Sick-Sicker model. The evolution of the cohort is simulated in one-year discrete-time cycles. Both costs and quality-adjusted life years (QALYs) are discounted at an annual rate of 3%.

State-transition diagram of the Sick-Sicker model. Healthy individuals can get Sick, die or stay healthy. Sick individuals can recover, transitioning back to healthy, can die, or stay sick. Once individuals are Sicker, they stay Sicker until they die.

Two alternative strategies exist for this hypothetical disease: a no-treatment and a treatment strategy. Under the treatment strategy, Sick and Sicker individuals receive treatment and continue doing so until they recover or die. The cost of the treatment is additional to the cost of being Sick or Sicker for one year. The treatment improves QoL for those individuals who are Sick but has no effect on the QoL of those who are sicker. To evaluate these two alternative strategies, we perform a CEA.

We assume that most of the parameters of the Sick-Sicker model and their uncertainty have been previously estimated and are known to the analyst. However, while we can identify those who are afflicted with the illness through obvious symptoms, we can not easily distinguish those in the Sick state from the those in the Sicker state. Thus, we can not directly estimate state-specific mortality hazard rate ratios, nor do we know the transition probability of progressing from Sick to Sicker. Therefore, we calibrate the model to different epidemiological data. We internally validated the calibrated model by comparing the predicted outputs from the model evaluated at the calibrated parameters against the calibration targets (Eddy et al. 2012, Goldhaber-Fiebert, Stout, and Goldie (2010)).

As part of the CEA, we conducted different deterministic sensitivity analysis (SA), including one-way and two-way SA, and tornado plots. To quantify the effect of parameter uncertainty on decision uncertainty, we conducted a probabilistic sensitivity analysis (PSA) and reported our uncertainty analysis results with a cost-effectiveness acceptability curve (CEAC), cost-effectiveness acceptability frontier (CEAF) and expected loss curves (ELC) (Alarid-Escudero et al. 2019). We also conducted a value of information (VOI) analysis to determine whether potential future research is needed to reduce parameter uncertainty. All steps of the CEA will be described using the different components of the framework.

Set up

This report is a supplementary material meant to guide you through the R code of a fully functional decision model to showcase the framework described by the Decision Analysis in R for Technologies in Health (DARTH) workgroup in the manuscript A need for change! A coding framework for improving transparency in decision modeling. The code of this analysis can be downloaded from GitHub (https://github.com/DARTH-git/Decision-Modeling-Framework). We recommend downloading the case-study files as a single .zip file containing all directories. Unzip the folder and save to your desired directory. The framework is divided into different directories, described in Table 1, that could be accessed from the RStudio project Decision-Modeling-Framework.Rproj. In this framework, you will find multiple directories as described in Table 1 of the main manuscript. We refer to the directory names of this framework and scripts stored in these directories using italic style. This report is created with Markdown and is located in the reports directory of the framework. The figures for the case-study can be found in the figs directory, data required to conduct some of the analyses of the different components are in the data directory and the R scripts with functions, are located in the functions directory. The main R scripts that conduct the analyses of the different components of the framework are stored in the R directory. In this document we do not show all the R code we refer to. Therefore, it is important to follow along while reading this document. To make sure you have all the required packages needed to run all the code, first install all the required packages by running the 00_general-setup.R script in the R directory.

source("R/app0_packages-setup.R")

01 Define model inputs

As described in the main manuscript, in this component we declare all model input variables and set their values. The R script running the analysis of this component is the 01_model-inputs.R file in the R directory.

The input to inform the values is divided in three categories: external, estimated, and calibrated. The majority of the Sick-Sicker model parameters are informed by external data. Only three parameter values need to be estimated using model calibration.

In this component, we start with the general setup of the model, specifying among others the time horizon, name and number of health states, proportion of the cohort in each of the different health states at the start of the simulation and discount rates. The next step is to specify the external parameters. The initial model parameter values and R variable names are presented in Table .

Description of the initial parameters with their R name and value of the Sick-Sicker model.

Parameter	R name	Value
Time horizon (\(n_t\))	n.t	75 years
Names of health states (\(n\))	v.n	H, S1, S2, D
Annual discount rate (costs/QALYs)	d.c/d.e	3%
Annual transition probabilities
- Disease onset (H to S1)	p.HS1	0.15
- Recovery (S1 to H)	p.S1H	0.5
- Disease progression (S1 to S2) in the time-homogenous model	p.S1S2	0.105
Annual mortality
- All-cause mortality (H to D)	p.HD	age-specific
- Hazard rate ratio of death in S1 vs H	hr.S1	3
- Hazard rate ratio of death in S2 vs H	hr.S2	3
Annual costs
- Healthy individuals	c.H	$2,000
- Sick individuals in S1	c.S1	$4,000
- Sick individuals in S2	c.S2	$15,000
- Dead individuals	c.D	$0
- Additional costs of sick individuals treated in S1 or S2	c.Trt	$12,000
Utility weights
- Healthy individuals	u.H	1.00
- Sick individuals in S1	u.S1	0.75
- Sick individuals in S2	u.S2	0.50
- Dead individuals	u.D	0.00
Intervention effect
- Utility for treated individuals in S1	u.Trt	0.95

Age-specific background mortality for healthy individuals is represented by the US population in 2015 and obtained from the Human Mortality database. This information is stored in the 01_all-cause-mortality.csv file in the data directory. Based on this .csv file a vector with mortality rates by age is created using the f.load_mort_data function in the 01_model-inputs_functions.R script. This function gives us the flexibility to easily import data from other countries or years.

print.function(f.load_mort_data) # print the function

## function (file = "data/01_all-cause-mortality.csv") 
## {
##     df.r.mort_by_age <- read.csv(file = file)
##     v.r.mort_by_age <- df.r.mort_by_age %>% dplyr::select(Total) %>% 
##         as.matrix()
##     return(v.r.mort_by_age)
## }

Another function in the 01_model-inputs_functions.R script, is the f.load_all_parms function. This function, which is actually using the f.load_mort_data function, loads all parameters for the decision model from multiple sources and creates a list that contains all parameters and their values.

print.function(f.load_all_params)  # print the function

## function (file.init = "data/01_init-params.csv", file.mort = "data/01_all-cause-mortality.csv") 
## {
##     df.params.init <- read.csv(file = file.init)
##     v.r.mort_by_age <- f.load_mort_data(file = file.mort)
##     l.params.all <- with(as.list(df.params.init), {
##         v.age.names <- n.age.init:(n.age.init + n.t - 1)
##         v.n <- c("H", "S1", "S2", "D")
##         n.states <- length(v.n)
##         v.s.init <- c(H = 1, S1 = 0, S2 = 0, D = 0)
##         l.params.all <- list(n.age.init = n.age.init, n.t = n.t, 
##             v.age.names = v.age.names, v.n = v.n, n.states = n.states, 
##             v.s.init = c(H = 1, S1 = 0, S2 = 0, D = 0), v.r.mort_by_age = v.r.mort_by_age)
##         return(l.params.all)
##     })
##     l.params.all <- c(l.params.all, df.params.init)
## }
## <bytecode: 0x7fa5c9797108>

The f.load_all_params function is informed by the arguments file.init and file.mort. The file.init argument is a string with the location and name of the file with initial set of parameters. The initial parameter values for our case-study are stored in the 01_init-params.csv file located in the data directory. The f.load_all_params function read this .csv file into the function environment as a dataframe called, df.params.init.

The file.mort argument is a string with the location and name of the file with mortality data. As described before, in our case-study this is the 01_all-cause-mortality.csv file. Within the f.load_all_parms function, the f.load_mort_data function is used to create a vector with mortality rates from the .csv data.

After loading all the information, the f.load_all_params generates a list called,l.params.all, including all parameters for the model including the general setup parameters and the vector of mortality rates. The function also stores the dataframe df.params.init with the initial set of parameters in the list. This is all executed in the in the 01_model-inputs.R script by running the code below.

l.params.all <- f.load_all_params(file.init = "data/01_init-params.csv",
                                  file.mort = "data/01_all-cause-mortality.csv")

For the Sick-Sicker model we do not have to estimated parameters, but we do have three parameters that need to be estimated via model calibration. In this stage of the framework, we simply set these parameters to valid “dummy” values that are compatible with the next phase of the analysis, model implementation, but are ultimately just placeholder values until we conduct the calibration phase. This means that these values will be replaced by the best-fitted calibrated values after we performed the calibration in component 3.

Using a function to create a list of base-case parameters to have all model parameters in a single object is very useful, because this object will have to be updated for the calibration and the different sensitivity analyses in components 3 and 5 of the framework, respectively. Below, we guide you through the components of the function.

02 Model implementation

In this component, we build the backbone of the decision analysis: the implementation of the model. This component is performed by the 02_simulation-model.R script. This file itself is not very large. It simply loads some packages, sources the input from component 01, sources the function f.decision_model that is used the capture the dynamic process of the Sick-Sicker example, runs this function and stores the output. The output of the model is the traditional cohort trace, describing how the cohort is distributed among the different health states over time, which is plotted at the end of this script. This trace will be used in many of the other components.

The function f.decision_model is defined in the 02_simulation-model_functions.R file. As described in the paper, constructing a model as a function at this stage facilitates subsequent stages of the model development and analysis, as these processes will all call the same model function, but pass different parameter values and/or calculate different final outcomes based on the model outputs. In the next part, we will describe the code within the function.

print.function(f.decision_model) # print the code of the function

## function (l.params.all, verbose = FALSE) 
## {
##     with(as.list(l.params.all), {
##         p.HDage <- 1 - exp(-v.r.mort_by_age[(n.age.init + 1) + 
##             0:(n.t - 1)])
##         p.S1Dage <- 1 - exp(-v.r.mort_by_age[(n.age.init + 1) + 
##             0:(n.t - 1)] * hr.S1)
##         p.S2Dage <- 1 - exp(-v.r.mort_by_age[(n.age.init + 1) + 
##             0:(n.t - 1)] * hr.S2)
##         a.P <- array(0, dim = c(n.states, n.states, n.t), dimnames = list(v.n, 
##             v.n, 0:(n.t - 1)))
##         a.P["H", "H", ] <- (1 - p.HDage) * (1 - p.HS1)
##         a.P["H", "S1", ] <- (1 - p.HDage) * p.HS1
##         a.P["H", "D", ] <- p.HDage
##         a.P["S1", "H", ] <- (1 - p.S1Dage) * p.S1H
##         a.P["S1", "S1", ] <- (1 - p.S1Dage) * (1 - (p.S1S2 + 
##             p.S1H))
##         a.P["S1", "S2", ] <- (1 - p.S1Dage) * p.S1S2
##         a.P["S1", "D", ] <- p.S1Dage
##         a.P["S2", "S2", ] <- 1 - p.S2Dage
##         a.P["S2", "D", ] <- p.S2Dage
##         a.P["D", "D", ] <- 1
##         m.indices.notvalid <- arrayInd(which(a.P < 0 | a.P > 
##             1), dim(a.P))
##         try(if (dim(m.indices.notvalid)[1] != 0) {
##             v.rows.notval <- rownames(a.P)[m.indices.notvalid[, 
##                 1]]
##             v.cols.notval <- colnames(a.P)[m.indices.notvalid[, 
##                 2]]
##             v.cycles.notval <- dimnames(a.P)[[3]][m.indices.notvalid[, 
##                 3]]
##             df.notvalid <- data.frame(`Transition probabilities not valid:` = matrix(paste0(paste(v.rows.notval, 
##                 v.cols.notval, sep = "->"), "; at cycle ", v.cycles.notval), 
##                 ncol = 1), check.names = FALSE)
##             if (verbose) {
##                 message("Not valid transition probabilities")
##                 stop(print(df.notvalid), call. = FALSE)
##             }
##         })
##         valid <- apply(a.P, 3, function(x) all.equal(sum(rowSums(x)), 
##             n.states))
##         if (!isTRUE(all.equal(as.numeric(sum(valid)), as.numeric(n.t)))) {
##             if (verbose) {
##                 stop("This is not a valid transition Matrix")
##             }
##         }
##         m.M <- matrix(0, nrow = (n.t + 1), ncol = n.states, dimnames = list(0:n.t, 
##             v.n))
##         m.M[1, ] <- v.s.init
##         for (t in 1:n.t) {
##             m.M[t + 1, ] <- m.M[t, ] %*% a.P[, , t]
##         }
##         return(list(a.P = a.P, m.M = m.M))
##     })
## }
## <bytecode: 0x7fa5cf368490>

The f.decision_model function is informed by the argument l.params.all. Via this argument we give the function a list with all parameters of the decision model. For the Sick-Sicker model, these parameters are stored in the list l.params.all, which we passed into the function as shown below.

l.out.stm <- f.decision_model(l.params.all = l.params.all) # run the function

This function itself has all the mathematical equations of the decision models coded inside. It starts by calculating the age-specific transition probabilities from all non-dead states based on the vector of age-specific mortality rates v.r.asr. These parameters will become vectors of length n.t, describing the probability to die for all ages from all non-dead states.

The next part of the function, creates an array that stores the age-specific transition probability matrices in each of the third dimension. The transition probability matrix is a core component of a state-transition cohort model (Iskandar 2018). This matrix contains the probabilities of transitioning from the current health state, indicated by the rows, to the other health states, specified by the columns. Since we have age-specific transition probabilities, the transition probability matrix is different at each cycle. These probabilities are only depending on the age of the cohort, and not on other events; therefore, we can generate all matrices at the start of the model. This results in n.t different age-specific matrices that are stored in an array, called a.P, of dimensions n.s x n.s x n.t. After initializing the array, it is filled with the transition probability stored in the list. When running the model, we can index the correct transition probability matrix corresponding with the current age of the cohort. We then added some sanity checks to make sure that the transition matrices and the transition probabilities are valid. The transition probability matrices stored in the array a.P, for the first three and last cycle, are shown below.

l.out.stm$a.P[, , 1:3] # show the first three time-points of a.P

## , , 0
## 
##            H        S1        S2           D
## H  0.8491385 0.1498480 0.0000000 0.001013486
## S1 0.4984813 0.3938002 0.1046811 0.003037378
## S2 0.0000000 0.0000000 0.9899112 0.010088764
## D  0.0000000 0.0000000 0.0000000 1.000000000
## 
## , , 1
## 
##            H        S1        S2            D
## H  0.8491513 0.1498502 0.0000000 0.0009985012
## S1 0.4985037 0.3938180 0.1046858 0.0029925135
## S2 0.0000000 0.0000000 0.9900597 0.0099402657
## D  0.0000000 0.0000000 0.0000000 1.0000000000
## 
## , , 2
## 
##            H        S1        S2           D
## H  0.8490910 0.1498396 0.0000000 0.001069428
## S1 0.4983976 0.3937341 0.1046635 0.003204853
## S2 0.0000000 0.0000000 0.9893570 0.010642959
## D  0.0000000 0.0000000 0.0000000 1.000000000

l.out.stm$a.P[, , l.params.all$n.t] # show it for the last cycle

##            H        S1         S2         D
## H  0.6055199 0.1068564 0.00000000 0.2876237
## S1 0.1807584 0.1427991 0.03795926 0.6384833
## S2 0.0000000 0.0000000 0.03365849 0.9663415
## D  0.0000000 0.0000000 0.00000000 1.0000000

By comparing these probability matrices, we observe an increase in the probabilities of transitioning to death from all health states.

After the array is filled, the cohort trace matrix, m.M, of dimensions n.t x n.s is initialized. This matrix will store the state occupation at each point in time. The first row of the matrix is informed by the initial state vector v.s.init. For the remaining points in time, we iteratively multiply the cohort trace with the age-specific transition probability matrix corresponding to the specific cycle obtained by indexing the array a.P appropriately. All the outputs and relevant elements of the decision model are stored in a list, called l.out.stm. This list contains the array of the transition probability matrix for all cycles t and the cohort trace m.M.

head(l.out.stm$m.M)    # show the top part of the cohort trace

##           H        S1         S2           D
## 0 1.0000000 0.0000000 0.00000000 0.000000000
## 1 0.8491385 0.1498480 0.00000000 0.001013486
## 2 0.7957468 0.1862564 0.01568695 0.002309774
## 3 0.7684912 0.1925699 0.03501425 0.003924648
## 4 0.7484793 0.1909659 0.05478971 0.005765035
## 5 0.7306193 0.1873106 0.07413838 0.007931783

tail(l.out.stm$m.M)    # show the bottom part of the cohort trace

##              H           S1           S2         D
## 70 0.009928317 0.0022433565 2.035951e-04 0.9876247
## 71 0.007153415 0.0015935925 1.311619e-04 0.9911218
## 72 0.005058845 0.0011174473 8.674540e-05 0.9937370
## 73 0.003460336 0.0007552206 5.436484e-05 0.9957301
## 74 0.002289594 0.0004937081 3.298632e-05 0.9971837
## 75 0.001475636 0.0003151589 1.985106e-05 0.9981894

Using the code below, we can graphically show the model dynamics by plotting the cohort trace. Figure shows the distribution of the cohort among the different health states at each time point.

Cohort trace of the Sick-Sicker cohort model

03 Model calibration

In this component, we calibrate unknown model parameters by matching model outputs to specified calibration targets. Specifically, we calibrate the Sick-Sicker model to match survival, prevalence and the proportion who are Sicker, among all those afflicted (Sick+Sicker). We used a Bayesian calibration approach using the incremental mixture importance sampling (IMIS) algorithm (Steele, Raftery, and Emond 2006), which has been used to calibrate health policy models (Raftery and Bao 2010, Menzies et al. (2017), Rutter et al. (2018)). Bayesian methods allow us to quantify the uncertainty in the calibrated parameters even in the presence of non-identifiability (Alarid-Escudero et al. 2018). This analysis is coded in the 03_calibration.R file. The target data is stored in the 03_calibration-targets.RData file. Similar to component 02, in the section 03.1 Load packages, data and functions we start by loading inputs and functions. In addition, we load the calibration targets data into the R workspace. In the next section, 03.2 Visualize targets, we plot each of the calibration targets with their confidence intervals.

In section 03.3 Run calibration algorithms, we set the parameters we like to calibrate to fixed values and test if the function f.calibration_out that produces model outputs corresponding to the calibration targets works. This function takes a vector of parameters that need to be calibrated and a list with all parameters of decision model and computes model outputs to be used for calibration routines.

print.function(f.calibration_out) # print the functions

## function (v.params.calib, l.params.all) 
## {
##     l.params.all <- f.update_param_list(l.params.all = l.params.all, 
##         params.updated = v.params.calib)
##     l.out.stm <- f.decision_model(l.params.all = l.params.all)
##     v.os <- 1 - l.out.stm$m.M[, "D"]
##     v.prev <- rowSums(l.out.stm$m.M[, c("S1", "S2")])/v.os
##     v.prop.S2 <- l.out.stm$m.M[, "S2"]/rowSums(l.out.stm$m.M[, 
##         c("S1", "S2")])
##     l.out <- list(Surv = v.os[c(11, 21, 31)], Prev = v.prev[c(11, 
##         21, 31)], PropSicker = v.prop.S2[c(11, 21, 31)])
##     return(l.out)
## }

This function is informed by two argument v.params.calib and l.params.all. The vector v.params.calib contains the values of the three parameters of interest and the list l.params.all contains all parameters of the decision model. The placeholder values are replaced by v.params.calib and with these values the model is evaluated. This is done by running the f.decision_model function, described in component 02. This results in a new list with output of the model corresponding to the parameter values in the v.params.calib. With this new decision model output, the overall survival, disease prevalence and the proportion of Sicker in the Sick and Sicker states are calculated. The estimated values for these epidemiological outcomes at different timepoints are combined in a list called l.out produced but the f.calibration_out.

Once we make sure this code works, we specify the calibration parameters in section 03.3.1 Specify calibration parameters. These include setting the seed for the random number generation, specifying the number of random samples to obtain from the calibrated posterior distribution, the name of the input parameters and the range of these parameters that will inform the prior distributions of the calibrated parameters, and the name of the calibration targets: Surv, Prev, PropSick.

In the next section, 03.3.2 Run IMIS algorithm, we calibrate the Sick-Sicker model with the IMIS algorithm. For this case-study, we assume a normal likelihood and uniform priors but other could be used. For a more detailed description of IMIS for Bayesian calibration, different likelihood functions and prior distributions, we refer the reader to the tutorial for Bayesian calibration by Menzies et al. (Menzies et al. 2017). We use the IMIS function from the IMIS package that calls the functions likelihood, sample.prior and prior, to draw samples from the posterior distribution (A. Raftery and Le Bao 2012). The functions are specified in the 03_calibration_functions.R file. For the IMIS function, we specify the incremental sample size at each iteration of IMIS, the desired posterior sample size at the resample stage, the maximum number of iterations in IMIS and the number of optimizers which could be 0. The function returns a list, which we call l.fit.imis, with the posterior samples, the diagnostic statistics at each IMIS iteration and the centers of Gaussian components (A. Raftery and Le Bao 2012). We store the posterior samples in the matrix m.calib.post.

We then explore these posterior distributions in section 03.4 Exploring posterior distribution. We start by estimating the posterior mean, median and 95% credible interval, the mode and the maximum-a-posteriori (MAP). To estimate the mode, we use the package modeest (Poncet 2018). All for these summary statistics are combined in a dataframe called df.posterior.summ. Table shows the summary statistics of the posterior distribution.

Summary statistics of the posterior distribution

	Mean	2.5%	50%	97.5%	Mode	MAP
p.S1S2	0.1077822	0.0966154	0.107619	0.1206907	0.1071694	0.1078278
hr.S1	2.7151807	1.1007592	2.608970	4.4139330	1.9654838	2.2970418
hr.S2	9.5435354	7.1631569	9.534398	11.8839979	9.2717728	9.8870386

In section 03.4.2 Visualization of posterior distribution, we generate a pairwise scatter plot of the calibrated parameters (Figure ) and a 3D scatter plot of the joint posterior distribution (Figure ). These figures are saved in the figs directory.

Joint posterior distribution

Pairwise posterior distribution of calibrated parameters

Finally, the posterior distribution and MAP estimate from the IMIS calibration are stored in the file 03_imis-output.RData. Storing this data as an .Rdata file allows to import the data in following sections without needing to re-run the calibration component.

04 Validation

In this section, we check the internal validity of our Sick-Sicker model before we move on to the analysis components. To internally validate the Sick-Sicker model, we compare the model-predicted output evaluated at posterior parameters against the calibration targets. This is all done in the 04_validation.R script by loading all previously described functions and the generated calibration targets.

In section 04.2 Compute model-predicted outputs, we compute the model-predicted outputs for each sample of posterior distribution as well as for the MAP estimate. We then use the function f.data_summary to summarize the model-predicted posterior outputs into different summary statistics.

print.function(f.data_summary)

## function (data, varname, groupnames) 
## {
##     require(plyr)
##     summary_func <- function(x, col) {
##         c(mean = mean(x[[col]], na.rm = TRUE), median = quantile(x[[col]], 
##             probs = 0.5, names = FALSE), sd = sd(x[[col]], na.rm = TRUE), 
##             lb = quantile(x[[col]], probs = 0.025, names = FALSE), 
##             ub = quantile(x[[col]], probs = 0.975, names = FALSE))
##     }
##     data_sum <- ddply(data, groupnames, .fun = summary_func, 
##         varname)
##     data_sum <- plyr::rename(data_sum, c(mean = varname))
##     return(data_sum)
## }
## <bytecode: 0x7fa5d2644170>

This function is informed by three arguments, data, varname and groupnames.

The computation of the model-predicted outputs using the MAP estimate is done by inserting the v.calib.post.map data into the previously described f.calibration_out function. This function creates a list including the estimated values for survival, prevalence and the proportion of sicker individuals at cycles 10, 20 and 30.

In sections 04.6 Internal validation: Model-predicted outputs vs. targets, we check the internal validation by plotting the model-predicted outputs against the calibration targets (Figure -). The generated plots are saved as .png files, which could be used in reports later on without the need of re-running the code.

Survival data: Model-predicted outputs vs targets

Prevalence data of sick individuals: Model-predicted output vs targets

Proportion who are Sicker, among all those afflicted (Sick + Sicker): Model-predicted output

05 Analysis

The analysis component is where the elements in components 1-4 are combined to answer the question(s) of interest given current information and to quantify the value of potential further research. Our framework separates the analysis in three subcomponents: 05a Deterministic analysis, 05b Uncertainty analysis and 05c Value of information analysis. For the Sick-Sicker case-study, we use all three subcomponents to conduct the CEA and to quantify the uncertainty of our decision. For procedures in the CEA, we rely on the R package dampack, which is available here: https://github.com/DARTH-git/dampack. Instructions for installing dampack are described in Appendix 0 provided in the app0_packages-setup.R script.

05a Deterministic analysis

In this subcomponent, we perform a deterministic CEA, followed by some deterministic sensitivity analysis, including one-way, two-way and tornado sensitivity analyses. The function script of this subcomponent, 05a_deterministic-analysis_function.R, contains the function f.calculate_ce_out. This function calculates costs and effects for a given vector of parameters using a simulation model. We need to run our simulation model using the calibrated parameter values, but the list we created in component 01 still contain the placeholder values for the calibrated parameters. This means we need to update these values by the calibrated values stored in the vector v.calib.post.map. The 00_general_functions.R script in the functions folder creates the f.update_param_list function that is able to update a list of parameters with new values for some specific parameters.

print.function(f.update_param_list)

## function (l.params.all, params.updated) 
## {
##     if (typeof(params.updated) != "list") {
##         params.updated <- split(unname(params.updated), names(params.updated))
##     }
##     l.params.all <- modifyList(l.params.all, params.updated)
##     return(l.params.all)
## }
## <bytecode: 0x7fa5c96cd6e0>

The first argument of the function, called l.params.all, is a list with all the parameters of decision model. The second argument, params.updated, is an object with parameters for which values need to be updated. The function returns the list l.params.all with updated values.

In the 05a_deterministic-analysis.R script we execute the f.update_param_list function for our case-study, resulting in the list l.params.basecase where the placeholder values for p.S1S2, hr.S1 and hr.S2 are replaced by the calibration estimates.

l.params.basecase <- f.update_param_list(l.params.all, v.calib.post.map) 

We use this new list as an argument in the f.calculate_ce_out function. In addition, we specify the willingness-to-pay (WTP) threshold value using the n.wtp argument of this function. This WTP value is used to compute a net monetary benefit (NMB) value. If the user does not specify the WTP, a default value of $100,000/QALY will be used by the function.

df.out.ce <- f.calculate_ce_out(l.params.all = l.params.basecase, 
                                n.wtp = 150000)
print.function(f.calculate_ce_out) # print the function

## function (l.params.all, n.wtp = 1e+05) 
## {
##     with(as.list(l.params.all), {
##         v.dwc <- 1/((1 + d.e)^(0:(n.t)))
##         v.dwe <- 1/((1 + d.c)^(0:(n.t)))
##         l.model.out.no_trt <- f.decision_model(l.params.all = l.params.all)
##         l.model.out.trt <- f.decision_model(l.params.all = l.params.all)
##         m.M_no_trt <- l.model.out.no_trt$m.M
##         m.M_trt <- l.model.out.trt$m.M
##         v.u_no_trt <- c(u.H, u.S1, u.S2, u.D)
##         v.u_trt <- c(u.H, u.Trt, u.S2, u.D)
##         v.c_no_trt <- c(c.H, c.S1, c.S2, c.D)
##         v.c_trt <- c(c.H, c.S1 + c.Trt, c.S2 + c.Trt, c.D)
##         v.tu_no_trt <- m.M_no_trt %*% v.u_no_trt
##         v.tu_trt <- m.M_trt %*% v.u_trt
##         v.tc_no_trt <- m.M_no_trt %*% v.c_no_trt
##         v.tc_trt <- m.M_trt %*% v.c_trt
##         tu.d_no_trt <- t(v.tu_no_trt) %*% v.dwe
##         tu.d_trt <- t(v.tu_trt) %*% v.dwe
##         tc.d_no_trt <- t(v.tc_no_trt) %*% v.dwc
##         tc.d_trt <- t(v.tc_trt) %*% v.dwc
##         v.tc.d <- c(tc.d_no_trt, tc.d_trt)
##         v.tu.d <- c(tu.d_no_trt, tu.d_trt)
##         v.nmb.d <- v.tu.d * n.wtp - v.tc.d
##         df.ce <- data.frame(Strategy = v.names.str, Cost = v.tc.d, 
##             Effect = v.tu.d, NMB = v.nmb.d)
##         return(df.ce)
##     })
## }
## <bytecode: 0x7fa5b5097c08>

After calculating the discount weights, this function runs the simulation model using the previously described function f.decision_model in the 02_simulatiomn-model_function.R script. Inside the function f.calculate_ce_out, the simulation model is run for both the treatment, l.model.out.trt, and no treatment, l.model.out.no_trt, strategies of the Sick-Sicker model. Running it for both treatment strategies is done for illustration purposes. In this case-study, the resulting cohort traces are identical and we could have executed it only once.

In the second part of the function we create multiple vectors for both the cost and effects of both strategies. These vectors multiply the cohort trace to compute the cycle-specific rewards. This results in vectors of total costs (v.tc) and total effects (v.tu) per cycle. By multiplying these vectors with the vectors with the discount weights for costs (v.dwc) and effects (v.dwe) we get the total discounted mean costs (tc.d_no_trt and tc.d_trt) and QALYs (tu.d_no_trt and tu.d_trt) for both strategies. These values are used in the calculation of the NMB. Finally, the total discounted costs, effectiveness and NMB are combined in the dataframe df.ce. The results for our case-study are shown below.

df.out.ce # print the dataframe 

##       Strategy     Cost   Effect     NMB
## 1 No Treatment 115239.8 20.02366 2888309
## 2    Treatment 214157.6 20.72299 2894291

This dataframe of CE results can be used as an argument in the calculate_icers function from the dampack package to calculate the incremental cost-effectiveness ratios (ICERs) and noting which strategies are weakly and strongly dominated. Table shows the result of the deterministic CEA.

df.cea.det <- calculate_icers(cost = df.out.ce$Cost, 
                              effect = df.out.ce$Effect, 
                              strategies = v.names.str)

Deterministic cost-effectiveness analysis results of the Sick-Sicker model comparing no treatment with treatment

Strategy	Cost	Effect	Inc_Cost	Inc_Effect	ICER
No Treatment	115239.8	20.02366	NA	NA	NA
Treatment	214157.6	20.72299	98917.83	0.6993333	141445.9

Finally, Figure shows the cost-effectiveness frontier of the CEA.

One-way sensitivity analysis results for the parameters c.Trt, p.HS1, u.S1 and u.Trt

We then conduct a series of deterministic sensitivity analysis. First, we conduct a one-way sensitivity analysis (OWSA) on the variables c.Trt, p.HS1, u.S1 and u.Trt and a two-way sensitivity analysis (TWSA) using the owsa_det and twsa_det functions in the 00_general_functions.R file. We use the output of these functions to produce different SA plots, such as OWSA tornado, one-way optimal strategy and TWSA plots (Figures ] - ).

One-way sensitivity analysis results for the parameters c.Trt, p.HS1, u.S1 and u.Trt

The optimal strategy with OWSA

The tornado plot

The two-way sensitivity results of the parameters u.Trt and u.S1.

05b Uncertainty analysis

In this subcomponent, we evaluate decision uncertainty by propagating the uncertainty through the CEA using probabilistic sensitivity analysis (PSA). Until now we used the parameter values as described in Table . However, we are uncertain about these values. Most of these input parameters are defined by probability distribution as described in Table .

Description of parameters with their R name and distribution.

Parameter	R name	Distribution
Annual transition probabilities
- Disease onset (H to S1)	p.HS1	beta(30, 170)
- Recovery (S1 to H)	p.S1H	beta(60, 60)
Annual costs
- Healthy individuals	c.H	gamma(shape = 100, scale = 20)
- Sick individuals in S1	c.S1	gamma(shape = 177.8, scale = 22.5)
- Sick individuals in S2	c.S2	gamma(shape = 225, scale = 66.7)
- Additional costs of sick individuals treated in S1 or S2	c.Trt	gamma(shape = 73.5, scale = 163.3)
Utility weights
- Healthy individuals	u.H	truncnorm(mean = 1, sd = 0.01, b = 1)
- Sick individuals in S1	u.S1	truncnorm(mean = 0.75, sd = 0.02, b = 1)
- Sick individuals in S2	u.S2	truncnorm(mean = 0.50, sd = 0.03, b = 1)
Intervention effect
- Utility for treated individuals in S1	u.Trt	truncnorm(mean = 0.95, sd = 0.02, b = 1)

In a PSA we sample the input parameter values from these distributions and we then run the model at each sample. In the file 05b_uncertainty-analysis_functions.R we created a single function, called f.generate_psa_params, that generates a PSA dataset for all the CEA input parameters. We specify the number of PSA samples via the n.sim argument. The function also accepts specifying a seed to allow reproducibility of the results.

print.function(f.generate_psa_params) # print the function 

## function (n.sim, seed = 20190220) 
## {
##     load("output/03_imis-output.RData")
##     n.sim <- nrow(m.calib.post)
##     set.seed <- seed
##     df.psa.params <- data.frame(m.calib.post, p.HS1 = rbeta(n.sim, 
##         30, 170), p.S1H = rbeta(n.sim, 60, 60), c.H = rgamma(n.sim, 
##         shape = 100, scale = 20), c.S1 = rgamma(n.sim, shape = 177.8, 
##         scale = 22.5), c.S2 = rgamma(n.sim, shape = 225, scale = 66.7), 
##         c.Trt = rgamma(n.sim, shape = 73.5, scale = 163.3), c.D = 0, 
##         u.H = rtruncnorm(n.sim, mean = 1, sd = 0.01, b = 1), 
##         u.S1 = rtruncnorm(n.sim, mean = 0.75, sd = 0.02, b = 1), 
##         u.S2 = rtruncnorm(n.sim, mean = 0.5, sd = 0.03, b = 1), 
##         u.D = 0, u.Trt = rtruncnorm(n.sim, mean = 0.95, sd = 0.02, 
##             b = 1))
##     return(df.psa.params)
## }
## <bytecode: 0x7fa5b4f7efe8>

The function returns the df.psa.input dataframe with a PSA dataset of the input parameters. With this dataframe we can run the PSA to produce distributions of costs, effectiveness and NMB. The PSA is performed by the 05b_probabilistic-analysis.R script. As shown in the code below, the df.psa.input dataframe is used by the f.update_param_list function to generate the corresponding list of parameters for the PSA. For each simulation, we perfrom three steps. First, the list of parameters is updated by the f.update_param_list function. Second, the model is executed by the f.calculate_ce_out function using the updated parameter list and third, the dataframes df.c and df.e store the estimated cost and effects, respectively. The final part of this loop is to satisfy the modeler when waiting on the results, by displaying the simulation progress.

for(i in 1:n.sim){ 
  l.psa.input <- f.update_param_list(l.params.all, df.psa.input[i,])
  df.out.temp <- f.calculate_ce_out(l.psa.input)
  df.c[i, ] <- df.out.temp$Cost
  df.e[i, ] <- df.out.temp$Effect
  # Display simulation progress
  if(i/(n.sim/10) == round(i/(n.sim/10),0)) {
    cat('\r', paste(i/n.sim * 100, "% done", sep = " "))
  }
}

We can plot the results using the plot function from dampack. Figure shows the CE scatter plot with the joint distribution of costs and effects for each strategy and their corresponding 95% confidence ellipse.

The cost-effectiveness plane graph showing the results of the probabilistic sensitivity analysis for the Sick-Sicker case-study.

Probabilistic cost-effectiveness analysis results of the Sick-Sicker model comparing no treatment with treatment.

Strategy	Cost	Effect	Inc_Cost	Inc_Effect	ICER
No_Treatment	115819.4	19.92611	NA	NA	NA
Treatment	214696.9	20.61257	98877.47	0.6864686	144037.9

Next, we perform a CEA using the previously used calculate_icers functions from dampack. Table shows the results of the probabilistic CEA. In addition, we plot a cost-effectiveness plane with the frontier, the cost-effectiveness acceptability curves (CEACs) and frontier (CEAF), expected Loss curves (ELCs) (Figure - ) (Alarid-Escudero et al. 2019). Followed by creating linear regression metamodeling sensitivity analysis graphs (Figure - ). All generated figures are shown below and stored to the figs folder (H. Jalal et al. 2013).

Cost-effectiveness frontier

Cost-effectiveness acceptability curves (CEACs) and frontier (CEAF)

Expected Loss Curves

One-way sensitivity analysis (OWSA)

Optimal strategy with OWSA

Tornado plot

Two-way sensitivity analysis (TWSA)

05c Value of information

In the VOI component, the results from the PSA generated in the probabilistic analysis subcomponent are used to determine whether further potential research is needed. We use the calc_evpi function from the dampack package to calculate the expected value of perfect information (EVPI). Figure shows the EVPI for the different WTP values.

evpi <- calc_evpi(wtp = v.wtp, psa = l.psa)

Expected value of perfect information

References

Alarid-Escudero, F, EA. Enns, KM. Kuntz, TL. Michaud, and H Jalal. 2019. “‘Time Traveling Is Just Too Dangerous’ But Some Methods Are Worth Revisiting: The Advantages of Expected Loss Curves Over Cost-Effectiveness Acceptability Curves and Frontier.” Value in Health In Press.

Alarid-Escudero, F, RF MacLehose, Y Peralta, KM Kuntz, and Enns EA. 2018. “Nonidentifiability in Model Calibration and Implications for Medical Decision Making.” Medical Decision Making 38 (7): 810–21. doi:10.1177/0272989X18792283.

Eddy, David M., William Hollingworth, J. Jaime Caro, Joel Tsevat, Kathryn M. McDonald, and John B. Wong. 2012. “Model transparency and validation: A report of the ISPOR-SMDM modeling good research practices task force-7.” Medical Decision Making 32 (5): 733–43. doi:10.1177/0272989X12454579.

Enns, E A, L E Cipriano, C T Simons, and C Y Kong. 2015. “Identifying Best-Fitting Inputs in Health-Economic Model Calibration: A Pareto Frontier Approach.” Medical Decision Making 35 (2): 170–82. doi:10.1177/0272989X14528382.

Goldhaber-Fiebert, JD, NK Stout, and SJ Goldie. 2010. “Empirically evaluating decision-analytic models.” Value in Health 13 (5): 667–74. doi:10.1111/j.1524-4733.2010.00698.x.

Iskandar, R. 2018. “A theoretical foundation for state-transition cohort models in health decision analysis.” PloS One 13 (12). Public Library of Science: e0205543–e0205543. doi:10.1371/journal.pone.0205543.

Jalal, Hawre, Bryan Dowd, François Sainfort, and Karen M. Kuntz. 2013. “Linear regression metamodeling as a tool to summarize and present simulation model results.” Medical Decision Making 33 (7): 880–90. doi:10.1177/0272989X13492014.

Menzies, Nicolas A., Djøra I. Soeteman, Ankur Pandya, and Jane J. Kim. 2017. “Bayesian Methods for Calibrating Health Policy Models: A Tutorial.” PharmacoEconomics 35 (6). Springer International Publishing: 613–24. doi:10.1007/s40273-017-0494-4.

Poncet, P. 2018. “Modeest: Mode Estimation.” https://CRAN.R-project.org/package=modeest.

Raftery, A, and L Bao. 2010. “Estimating and Projecting Trends in HIV/AIDS Generalized Epidemics Using Incremental Mixture Importance Sampling.” Biometrics 66 (4): 1162–73.

Raftery, Adrian, and Le Bao. 2012. IMIS: Increamental Mixture Importance Sampling. https://CRAN.R-project.org/package=IMIS.

Rutter, C, J Ozik, M DeYoreo, and Collier N. 2018. “Microsimulation Model Calibration using Incremental Mixture Approximate Bayesian Computation.” arXiv, no. april: 1–20. https://arxiv.org/abs/1804.02090v3.

Steele, R, A Raftery, and M Emond. 2006. “Computing Normalizing Constants for Finite Mixture Models via Incremental Mixture Importance Sampling (IMIS).” Journal ofComputational and Graphical Statistics 15 (3): 712–34. doi:10.1198/106186006X132358.